Spinon-holon interactions in an anisotropic t-J chain: a comprehensive study 



O 

o 

(N 

c/3 



i 

u 
+-> 

CO 

-4-* 



i 

C 

O 

o 



> 
in 
in 

cn 

in 
o 
i> 
o 



Jurij Smakov, A. L. Chernyshev, Steven R. White 
Department of Physics and Astronomy, University of California Irvine, Irvine, California 92697, USA 

(Dated: February 1, 2008) 

We consider a generalization of the one-dimensional t-J model with anisotropic spin-spin interac- 
tions. We show that the anisotropy leads to an effective attractive interaction between the spinon 
and holon excitations, resulting in a localized bound state. Detailed quantitative analytic predic- 
tions for the dependence of the binding energy on the anisotropy are presented, and verified by 
precise numerical simulations. The binding energy is found to interpolate smoothly between a finite 
value in the t-J z limit and zero in the isotropic limit, going to zero exponentially in the vicinity 
of the latter. We identify changes in spinon dispersion as the primary factor for this non-trivial 
behavior. 

PACS numbers: 71.10.Fd, 71.10.Li, 75.10.Pq, 75.40.Mg 



I. INTRODUCTION 

One-dimensional (ID) lattice models of strongly cor- 
related fermions and bosons have traditionally been an 
object of intense theoretical studies. The reason for such 
interest is twofold. First, such models are relevant for the 
description of many real physical systems, such as mate- 
rials with strong uniaxial anisotropy, optical lattices, and 
quantum nanowires. Second, there is a number of the- 
oretical methods, unique to one dimension, which allow 
either exact (solution via Bethe Ansatz) or quasi-exact 
(bosonization, various renormalization group schemes) 
treatment of the models in question^ 

Out of the vast variety of ID models of strongly cor- 
related fermions, the one known as the t-J model clearly 
stands out as simple, yet remarkably versatile. It cap- 
tures both the ability of the particles to hop from one 
site to another, and the spin-spin interactions between 
them. By tuning the ratio of the coupling constants and 
the doping level, it may be used to describe many ID 
systems, ranging from non-interacting mobile fermions 
to Hcisenberg spin chains. Furthermore, it also repre- 
sents a physically relevant limit of another ID model of 
paramount importance - the Hubbard model. 

In the one-dimensional t-J model spin and charge dy- 
namics are independent, leading to the well-known effect 
of spin-charge separation^ the splitting of the electron 
(hole) into spinon and holon elementary excitations that 
carry only spin and only charge, respectively. This may 
be observed already at the single-hole doping level. In 
that case the low-energy spectrum of the t-J model has 
been extensively studied in the pasti^^i Recently it 
has been also shown that spinon and holon excitations are 
affected by effective attractive interaction which, how- 
ever, does not result in their binding or pairing. 6 This 
seeming controversy encouraged us to explore in detail 
the nature of spinon-holon interactions. In this paper 
we consider the t-J model as a limiting case of a more 
general model, which has anisotropic (XXZ-tike) spin- 
spin interactions. In this model the effective attractive 
spinon-holon interactions naturally emerge, leading to a 
spinon-holon bound state. We present detailed quantita- 



tive analytical predictions for the behavior of the binding 
energy as a function of anisotropy and the implications 
of this physical picture for the isotropic case, and verify 
them with precise numerical simulations, using exact di- 
agonalization (ED) and density-matrix renormalization 
group (DMRG) on systems of up to 23 and 128 sites, 
respectively. An experimental test of our work could 
come from photemission studies in insulating spin-chain 
materials with the Ising anisotropy, such as CSC0CI3, 
CsCoBr3, and others.— We note that in the past such ex- 
periments in the isotropic Heisenberg spin-chain material 
of the cuprate family, SrCuOs, have provided direct ex- 
perimental evidence of spin-charge separation in the real 
t-J model-like system^ 

The one-dimensional t-J model is defined by a Hamil- 
tonian H = ^ with 
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where c a i annihilates a fermion with spin a on site i, 
rii is the fermion number operator on site i, and Sj is 
the fermion spin operator. Periodic boundary conditions 
(BCs) are assumed. The Hilbert space, in which the 
Hamiltonian |1| acts, is restricted to a subspace without 
any doubly-occupied sites. At half-filling (one fermion 
per site) no particle hopping is possible, so the model is 
reduced to an isotropic Heisenberg model of interacting 
spins, with an antiferromagnetic (AF) ground state (GS). 
Doping it, even with a single hole, leads to spin-charge 
separation, which is manifested by the splitting of quasi- 
particle peaks in the excitation spectrum into two differ- 
ent sets, with energies scaling with t or J, respectively^ 
In order to study the spinon-holon interaction as a 
function of anisotropy, we consider a generalization of 
the t-J model ([!]) with a Hamiltonian Hij of the form 
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Here ■ Sf = SfSf + SfSV, and parameter a con- 
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trols the anisotropy of spin-spin interactions. The origi- 
nal isotropic t-J model is recovered by setting a — 1. 

The a = limit of Hamiltonian (JSJ) is known as t-J z 
model. Its GS in the undoped state is an Ising antifer- 
romagnet, and the effect of doping it with a single hole 
is easy to understand (see Fig. [T] for an illustration) . It 
results in a creation of a spinon-holon bound state due 
to an effective attraction between the immobile spinon 
(the Hamiltonian does not contain any spin-flipping term 
which would allow it to propagate) and a free holon^ 
The binding energy can then be calculated analytically: 



A = 2t 1 - y/1 + (J z /4ty 



(3) 
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We present several different methods to obtain this result 
in the Appendix. Setting a to a non-zero value presents 
three distinct possibilities. First of all, it is possible that 
any non-zero value of a immediately destroys the bound 
state, so a = is the only singular point in the phase di- 
agram with a finite A. Second, there is a possibility that 
A varies smoothly with a, interpolating between the fi- 
nite value at a = and zero value in the isotropic case. 
Finally, the A (a) dependence can go to zero at some non- 
trivial critical value < a c < 1. Out of these possibilities 
the first one appears to be the least likely one, as it is 
intuitively clear that small transverse spin-spin interac- 
tion cannot immediately destroy the bound state. While 
at a ^ the spinon will become mobile, for small a it 
is still going to be too "massive" , compared to virtually 
free holon. We cannot unequivocally rule out the last 
option (binding becomes too weak to be detected numer- 
ically near the isotropic limit), but we argue that in this 
regime the spin background is Ising-like, with long-range 
spin order for any anisotropy a < 1. This fact strongly 
suggests that the only anisotropy-driven critical point in 
the system is at a = 1. To confirm this hypothesis and 
carefully examine the remaining option, a detailed in- 
vestigation of the binding energy as a function of a is 
required. Such an investigation is the main topic of this 
paper. 

We have chosen the representative value of J z /t = 4.0 
for most of our calculations, after confirming that the 
results at other values of I < J z /t < 8 are qualitatively 
similar. We also present the final results for the binding 
energy for J z /t = 1.0. In general, we do not expect any 
qualitative difference for any other J z /t value as Eq. © 
gives A < for any J z . The choice of J z /t — 4.0 was 
made mostly to optimize the numerical accessibility of 
the binding energy in a wider range of a. 

The rest of the paper is organized as follows. In sec- 



FIG. 2: Lowest ED energies E s (k) - E a (0) vs k in L = 21 
chain with zero holes at J z /t — 4.0 (spinon dispersion). Solid 
lines show spinon dispersion for an infinite system obtained 
from BA.— All the energies are in units of t. 



tion |TT] we present our numerical results, discussing in 
detail the finite-size effects of the data and the proce- 
dure for extrapolation to the infinite system size. Section 
IIIII contains the theory for the binding energy, based on 
Bethe-Salpeter equation. We summarize our results in 
section HVl and present three different ways to derive the 
expression for the binding energy of t-J z model ([3]) in the 
Appendix. 



II. NUMERICAL RESULTS 

We have used the ED and DMRG techniques to cal- 
culate the ground state energies (GSEs) of the model for 
different system sizes and doping levels. This informa- 
tion was then used to extract the binding energy of a 
spinon-holon state in the infinite size limit. 

In ED we start by considering a subset of states of a 
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FIG. 1: (Color online). A hole in the Ising AF background 
(circle), moved by four sites from origin. The location of 
immobile spinon is indicated by the dashed box. 



FIG. 3: Lowest ED energies E h (k) - E h (0) vs k in L = 21 
chain with one hole at J z /t — 4.0 (holon dispersion). Dashed 
lines are guides to the eye. All the energies are in units of t. 
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system of size L with given total hole number n and to- 
tal S z . To take advantage of the translational symmetry, 
we then use these to construct a basis out of eigenstates 
of the translation operator with a given momentum k. 
Finally, the Hamiltonian matrix in this reduced basis is 
constructed, and its lowest eigenvalue is calculated itera- 
tively using the Lanczos algorithm. The implementation 
of every step in the procedure is described in detail in 
Ref. Ojl- With these techniques we were able to calcu- 
late the GSEs of systems of up to 23 sites with ED. Us- 
ing DMR G 13 ' 14 we have calculated GSEs of systems of 
up to L — 128 sites using periodic boundary conditions 
(PBCs), which greatly increases the numerical effort re- 
quired. Up to m = 1400 states per block were kept in 
the finite system method, with corrections applied to the 
density matrix to accelerate convergence with PBCs. 15 

We have carefully tested our algorithms by compar- 
ing the results of ED and DMRG for different system 
sizes, both with and without the hole. We have also 
compared the ED energies with the independent results 
for the GSEs of the XX Z models In all cases agreement 
to at least 7 decimal places was achieved. 

The elementary excitations of the model may be stud- 
ied by looking at systems of different sizes with either 
no or one hole. In the case of an odd number of sites 
and no holes the PBCs are frustrating, corresponding to 
creation of a frustrated ferromagnetic link - a spinon ex- 
citation. The lowest energies for each momentum sector 
for a chain of 21 sites with PBC at different anisotropics 
are shown in Fig. [21 This gives us the spinon disper- 
sion, which evolves from completely flat in the Ising case 
a = to quasi-relativistic in the isotropic case a = 1. 
Solid lines in Fig. [2] show the exact BA result for the 
spinon spectrum in the XX Z-mo&el: 11 
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k 2 sin 2 q. 



(4) 



Here c/J z = K\J\ — a 2 /ir, and k is determined from 
the condition ttK' /K — cosh -1 (1/a), where K = K(k) 
and K' = K(y/l — n 2 ) are complete elliptic integrals of 
the first kind. The lowest spinon energy is attained at 
q = ±ir/2 for any a > 0. 

A configuration with odd number of sites and one hole 
contains a "pure" holon, which can propagate through 
the system without disturbing the otherwise perfect AF 
background. Typical holon dispersions, obtained by mea- 
suring the energies of a 21-site chain with one hole, are 
presented in Fig. [3J The holon's minimum energy de- 
pendence on a is non-trivial. At a — the system has 
unique lowest energy point at momentum zero. This is 
only true for a finite system though, as in the infinite 
system this energy would be degenerate with the one at 
k = 7T. However, since we do not have a reciprocal space 
point exactly at k = ir, for a finite system the energy 
at the momentum points closest to k — tt is somewhat 
higher. This mismatch is an important source of finite- 
size effects in our measurement, as we discuss below. As 
the anisotropy a is increased, the energy at k — ir de- 



0.2 
0.4 
^-0.6 
0.8 
-1 



-1.2 r 

-1.4 r 




T 



A 



o~e~- e— e— e— o--gt=# 
: d i 



7l/4 



i 

nil 



o-o a=0.0 
Q^ia=0.1 
a-a a =0.5 
v-^a=1.0 



371/4 



71 

k 



5ti/4 3tc/2 77t/4 



2% 



FIG. 4: Lowest ED energies E p (k) — E p (0) vs k in L = 20 chain 
with one hole at J z /t — 4.0 (spinon-holon pair dispersion). 
Dashed lines are guides to the eye. 



creases, so the ground state switches from the k — to 
k = 7r sector at some finite intermediate value. Notably, 
these observations are in stark contrast with assumptions 
by Shiba and Ogata,— who claim that the k = and 
k = 7r energies are going to be degenerate for any finite 
system in the isotropic case (they are degenerate in an in- 
finite system though). Furthermore, their interpretation 
of the holon dispersion (presented in Fig. 6 of Ref. y) is 
somewhat misleading: they attribute the double-peaked 
structure of the dispersion to "strong antiferromagnetic 
correlations" . Our calculations confirm that the holon 
dispersion close to the k = and k = tt points may be 
very well fitted with a simple cosine dispersion of a free 
particle. This indicates that the characteristic double- 
peaked shape is formed by two different holon branches, 
centered at k — and k — tt. As one moves away from 
these points towards k = ir/2, the energy of the excita- 
tions grows, eventually making the creation of a spinon- 
antispinon pair energetically favorable, as suggested in 
Ref. |6|. That results in mixing of the two holon branches, 
which leads to a rounding of the dispersion peaks. 

Finally, an even-sized system with one hole corre- 
sponds to a situation where both spinon and holon are 
present. The GSE as a function of k for the system con- 
taining a spinon-holon pair is presented in Fig. 2] 

We can measure the energies of an interacting spinon- 
holon pair, as well as those of individual spinon and holon 
excitations, by taking the GSE of a corresponding config- 
uration and subtracting the extensive part of the energy 
e a L, where e a is the energy per site of an infinite XX Z 
chain, known from Bethe AnsatzJi That way we can ob- 
tain the spinon, holon, and spinon-holon pair energies for 
a set of different system sizes. After extrapolating to the 
infinite system size, the corresponding energies E Sl E/- t , 
and E p can be used to calculate the binding energy of 
the spinon-holon state in an infinite system as 

A = E p - E s - E h . (5) 

We will refer to this approach as "method A" . 
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In order to obtain an accurate estimate for the exci- 
tation energies in the infinite size limit, we have to deal 
with a variety of finite-size effects. The lifting of de- 
generacy in holon dispersion mentioned above is one of 
them. It turns out that its effect on the resulting GSE 
depends on whether the system size L (or L — 1 if L is 
odd) is divisible by 4 or not, so we will refer to these two 
data branches as 4-even and 4-odd, respectively. Such 
a mod (4) dependence has been extensively discussed in 
the literature (see Ref. and references therein). In the 
4-even branch the energy ek=o provides an upper bound 
for the true GSE, while efc =7r serves as a lower bound, 
and the bounds are reversed for the 4-odd case. Another 
source of finite-size corrections is the incommensurabil- 
ity of the momentum space points in the systems of odd 
size. For example, it can be seen from Fig. 03 that for 
spinons the GS corresponds to momentum tt/2 in the 
L = oo limit. However, for any finite-sized system with 
odd L there will be no reciprocal space point fc = tt/2, 
instead the GSE will occur at one of the nearest points 
with momentum fc — it/2±Sl, where with Sl = tt/L. As 
the system size is increased, Sl will go to zero, and the 
GSE will drift towards its infinite-!/ limiting value. Sim- 
ilarly, we cannot directly measure the GSE for holons 
(Fig. 02) at k = 7r. All these factors lead to a highly 
non-trivial finite size dependence. As an example, Fig. 
[5] shows the size dependence of the raw holon energies. 
To get a meaningful extrapolation the separate analysis 
of 4-even and 4-odd branches, which contain only half of 
the original points, is required. 

The situation with incommensurate fc-space points can 
be improved by imposing twisted boundary conditions on 
the models A boundary twist leads to the translation 
of the points in fc-space, but does not affect the energy 
spectrum. Thus, by adjusting the twist one can shift 
the fc-point with anticipated minimum energy from an 
incommensurate location in fc-space to an accessible one. 
In case of holons such a procedure is particularly simple, 
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FIG. 5: Raw holon GS energies for J z /t = 4.0 and differ- 
ent anisotropics a as a function of inverse system size 1/L. 
Dashed lines are guides to the eye. 
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FIG. 6: Ground state energy as the function of phase (j) f° r 
L — 16, J z /t — 4.0, a = 0.5. Dashed line is a quadratic fit. 
Inset shows the comparison of raw staggered pair energy data 
(open circles) and phase-corrected data (solid circles) for the 
same J z and a, and different system sizes. 



since we have to use a phase shift of 7r to move the fc = tt 
point of the original model to momentum fc = for a 
model with boundary twist. Such a phase shift is readily 
implemented just by switching the sign of the hopping 
constant t to the opposite one. That enabled us to re- 
construct the points, lost due to the splitting into 4-even 
and 4-odd branches, by complementing the holon GSE 
data with measurements performed on the model with 
t = — I , as shown in Table HI 

The spinon-holon pair GSE data also suffer from the 
fc-mismatch, as the GS is achieved at an incommensurate 
fc-point. 19 In principle, the same procedure may be ap- 
plied to improve the pair energy data. There, however, 
the phase shift <j) needed to shift the energy minimum 
to an accessible momentum point is size-dependent, so it 
has to be determined individually for every data point. 
By replacing t in (03) by te 1 ^ and tuning the phase shift <j>, 
we were able to measure the total energy of the system 
as a function of <p using ED. An example of such depen- 
dence is presented in Fig. 03 One remarkable feature of 
this dependence is that it is very well fit by a quadratic 
polynomial, so it is sufficient to know the energy at two 
different non-zero values of <j> to recover the "true" low- 
est energy at the minimum with excellent accuracy. The 
inset of Fig. 03 shows dramatic improvement of the data 



Branch 



4-even 



4-odd 



mod(L - 1,4) 



sign(t) 



-1 



TABLE I: Splitting of holon energy data for different sizes L 
and different signs of the hopping constant t into the 4-even 
and 4-odd branches. 
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due to the phase-induced correction. While such bound- 
ary conditions can be readily handled by ED, our DMRG 
code required extensive modifications to support them. 
Thus, in this work we perform the extrapolations using 
only the raw spinon-holon pair GSE data, split into 4- 
even and 4-odd branches. Further improvement of the 
precision of our results by using phase-adjusted data is 
possible. 

After the data for holon and pair are split into such 
branches, we need to extrapolate them to the L = oo 
limit, using a reasonable fitting form. From the holon 
(Fig. [7|) and pair (Fig. \§§ excitation energy data it 
is evident, that it has a complicated size dependence, 
which cannot be adequately described by a polynomial. 
Clearly, at large L the difference between the limiting 
value and the data points drops exponentially with in- 
creasing L. Incidentally, the size dependence for the 
GSE of the XXZ model, deduced by Woynarovich and 
de Vega (WdV) from BA^ is also dominated by an ex- 
ponential factor exp(— XL). That inspired us to attempt 
fitting the holon and pair data with the functional form 



El — E n 



-XL 



Pn(l/L), 



(6) 



where P n (x) is a polynomial of order n (n < 4) in x with 
adjustable coefficients. Not only does it work remarkably 
well for both holons and pairs (extrapolations are shown 
in Fig. [7] and Fig. [5] with dashed lines), but the values of 
the coefficient A we have found by keeping this parameter 
free in our holon fits provide an excellent match to the 
analytic values found by WdV for the XXZ model (their 
comparison is presented in Fig. [5]). Therefore, we have 
assumed that for holons the A values found by WdV are 
either exact, or a very good approximation. Thus, we 
used them in our holon fits, reducing the total number 
of free parameters by one. For pairs the values of the A 
parameter did not correlate with WdV results at all, so 
it had to be kept as a free parameter in the fit. 

When the parameter A is sufficiently large (for a < 
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FIG. 8: Comparison of the analytic values of A from Ref. yj| 
(dashed line) with the values obtained by fitting the holon 
data with form ((6]) and keeping A a free fitting parameter 
(symbols). 



0.3), the exponential factor in © makes the asymptotic 
approach to the infinite value very rapid, allowing us to 
simply adopt the energy value for the largest available 
size as the infinite-size limiting value. Increasing a re- 
sults in decreasing A, which pushes the onset of the expo- 
nential size dependence to larger and larger system sizes. 
In this regime the extrapolation using form ^ must be 
used. Around a ~ 0.5 parameter A becomes compara- 
ble with the inverse of the maximum available system 
size. At higher anisotropies the onset of the exponential 
behavior in size dependence takes place at characteristic 
sizes, not accessible by our calculations (as can be seen 
on the lower right panel of Fig. [5]), making the precise 
extrapolation of the pair excitation energy impossible. 
We can improve the accuracy of the extracted infinite 
size value -E^ by noting that both for holons and pairs 
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FIG. 7: The 4-even (circles) and 4-odd (squares) branches of 
holon energy data for J z /t = 4.0 and different anisotropies a. 
The "good" branch is used to extrapolate to L — oo, using 
form (O with n = 4 and the WdV values of A from Ref. [H 



FIG. 9: The 4-even (circles) and 4-odd (squares) branches of 
spinon-holon pair energy data for J z /t = 4.0 and different 
anisotropies a. The "good" branch is used to extrapolate to 
L — oo, using form (0 with n = 4. 
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Branch 


mod(L, 4) 


sign(t) for E\_i 


sign(i) for E l L+1 


Bl 





-1 


+1 


B2 





+1 


-1 


S3 


2 


-1 


+1 


54 
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+1 


-1 



TABLE II: Subdivision of the numerical data for A(L) into 
different branches due to finite size effects in method B. 



one of the branches is always more "well-behaved" than 
the other one. For example, non-uniform behavior of the 
4-even branch for holons can be seen on the lower left 
panel of Fig. [7] (it peaks slightly around 1/L = 0.025), 
and on the upper panels of Fig. [9] for the 4-odd pair 
branch. This non-uniformity of the "bad" branch usu- 
ally results from the GS switching from one momentum 
sector to a different one as a function of L. In our anal- 
ysis we have used only the extrapolations obtained with 
the "well-behaved" branch - 4-odd for holons and 4-even 
for pairs. The energy for the spinon excitations can, in 
principle, be extracted from the numerical data in a simi- 
lar way. However, to further improve our results, we have 
used the analytic expression for the spinon excitation en- 
ergy E s = w g=7r /2, given by Eq. ([4J, thus eliminating the 
finite-size effects from the spinon GSE completely 

Finally, the binding energy results for L — oo obtained 
by method A using ([5]) for J z /t — 4.0 and J z /t = 1.0 are 
presented in Fig. [TUJ These data are of high-precision 
for a < 0.5. We estimate the maximum relative error 
of the resulting binding energy by studying the quality 
of the fits and the variation of Aoo depending on the fit 
type. At a = 0.5 the error does not exceed 3% (10%) 
for J z /t = 4.0 (J z /t — 1.0) and becomes negligible very 
rapidly for smaller values of a. For a — 0.6 the error is of 
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FIG. 10: (Color online). Binding energy A as a function of 
a for Jz/t = 4.0 and J z /t — 1.0 (inset). Data includes the 
theoretical prediction (solid line), numerical results from ED 
and DMRG data obtained by method A (circles) and method 
B (diamonds). Dashed line shows the linear approximation 
(1251) . valid at small a. 




1/L 1/L 

FIG. 11: Size dependence for different branches of binding en- 
ergy A(L) at J z /t = 4.0. For branch definition see Table [Til 
Solid lines are guides to the eye, dashed curve is the polyno- 
mial extrapolation of the "good" branch B2. Legend applies 
to all panels. 

the order of 10% (100%) for J z /t = 4.0 {J z /t = 1.0), and 
for a — 0.7 it exceeds 100% for both representative values 
of J z . As mentioned before, we do not consider the bind- 
ing energy results for a > 0.6 to be reliable due to issues 
with pair energy extrapolation. Also, at larger values of 
a the value of the binding energy becomes comparable 
with the accuracy of our DMRG method (about 10" 7 to 
10~ 8 absolute precision, leading to about 10~ 4 t accuracy 
of the binding energy at L = 128 and J z /t = 4.0), im- 
posing a natural limitation on the quality of the data. 

An alternative way to extrapolate the binding energy 
to the infinite size limit is to calculate it for every sys- 
tem size L individually, and then do the extrapolation of 
the resulting size dependence to L = oo. In this method 
(referred to as "method B" ) we intentionally avoid using 
any BA results, to see whether the reliable binding en- 
ergy data may be obtained based on the numerical results 
alone. One could hope that the finite-size effects of vari- 
ous components entering the binding energy may cancel 
out, allowing the extrapolation to the infinite-size limit 
using a simple polynomial in 1/L, instead of an expo- 
nential. This approach, not depending on the theoretical 
results, provides an important validity test for the results 
of method A. 

Since holon and spinon GSEs are only available for 
odd L, and the spinon-holon pair ones only for even L, 
we define the finite size binding energy for an even size 
L as 

A(L)=El + El-[E° L _ 1 +El +1 +El_ 1 +El +1 }/2, (7) 

where E\ is the ground state of a system with L sites, 
doped with h holes. This expression is analogous to ([5]): 
sum of first two terms corresponds to the pair energy, 
while + E° L+1 )/2 and {E 1 L _ 1 + El +1 )/2 represent 

the average energy of a system of size L with a spinon 
and holon, respectively. Again, due to staggering of the 
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GSEs, binding energy data splits into a 4-even and 4- 
odd branches, depending on whether L is divisible by 4 or 
not. However, in this case we have an additional freedom 
of choosing the sign of t for holon energies E\_ x and 
E\j rl . Taking this into account results in 4 different data 
branches, defined in Tabic HT1 The remaining possibilities 
of using the same sign of t both for E\_ x and E\ +1 have 
been discarded as obviously suboptimal. 

The size dependence of different data branches for dif- 
ferent anisotropies is presented in Fig. [TTJ Due to split- 
ting, the number of points in each branch is pretty small, 
so we have used a polynomial of maximum possible de- 
gree (one less than the number of points) to perform the 
extrapolation to the infinite system size. From our previ- 
ous experience we know that the "good" holon branch in 
method A corresponds to branch B2, therefore we used 
the extrapolated value from this branch as our final re- 
sult for the binding energy. As can be seen in Fig. [TT)I 
the data from methods A and B are in excellent agree- 
ment in the range of a, where its calculation is reliable. 
However, we have found that the method B data always 
has a larger relative error than method A, mainly due to 
the size-dependence of the spinon component, eliminated 
in method A. 



III. THEORETICAL RESULTS 

Because in the Ising limit the spinon is impurity-like, 
the spinon-holon binding at a = can be solved in a 
number of ways (see Appendix) to give Eq. ([3]). For 
finite a the binding energy of the spinon-holon state may 
be calculated analytically by finding the poles of the two- 
particle scattering amplitude T. We will use shorthand 
notation q = (g, lo) and k = (fc, e) to denote the momenta 
and energies of spinon and holon, respectively. 

Scattering amplitude T obeys the Bethe-Salpeter 
equation,— presented in diagrammatic form in Fig. 1121 
Generally, it depends on both the incoming q, k and out- 
going q, k 2-momenta of spinon and holon. Using that 
momentum and energy are conserved, k + q = k + q = 
P = (P, E) , we may write it as 

r P (q,q) = V< M -+ / U 9Y G q ,G P _ q ,r P (q',q). (8) 

■V 

Here V q , q i is the spinon-holon interaction, and G h ^ is the 
holon (spinon) Green's function. A shorthand notation 





FIG. 12: Bethe-Salpeter equation for the spinon-holon scat- 
tering amplitude (circle). Spinons (holons) are shown by 
dashed (solid) lines. 



f = r°° V with V . = C *f is used. In 
the vicinity of the pole of T, Eq. ([5]) should reduce to a 
homogeneous integral equation with T whose dependence 
on one of the momenta q is only parametric and can be 
dropped^ 



Tp(q) = / F^G q ,G P _ q ,r P (q') 

J a' 



(9) 



The holon and spinon create a bound state if this integral 
equation has a solution. By introducing a function 



Xp(<?) = / G q G£_ q r P (q), 

J UJ 



(10) 



multiplying both sides of (O by G q G P _ q , and integrating 
over lo, we arrive at 



xp(q) 



°q U P q 



Evaluation of the first integral on the rhs requires knowl- 
edge of the spinon and holon Green's functions. For now 
we will just assume that they are free particles with some 
dispersions to q and e^, the specific form of which is to be 
determined: 



C s = 



iS 1 



G 



t - e k + i8 



With this assumption the integral is trivially done, yield- 
ing the final form of the Bethe-Salpeter equation for 

Xp(«): 



Xp(q) 



1 



E - ep- q - UJ, 



E^xp(?') 



(12) 



i' 



From this equation it is clear that xp(q) is nothing 
but the pair wavefunction and the equation (|12p is the 
Schrodinger equation for it in integral form. The pair 
energy E may be thought of as the binding energy A, 
measured relative to the lowest energies of the particles 
£o = min[tk\ and ujq = min[u>q\: 



E = A + e +uj . 



(13) 



In the Ising limit e k — ~2tcosk, eg = —2t, uo q = loq = 
J 2 /2, and Vq. q > = —uj , so Eq. (TT2")) is readily solved by 



Xp(<?) = 



C 



E - e P - 



(14) 



yielding a dispersionless (P-independent) bound state 
with A given by ([3]). From general considerations, the 
binding energy in ID should scale as —V 2 m, where V is 
interaction strength and m is the particle mass. In the 
Ising case this gives A^ —J^/t, in agreement with the 
exact result ©. 
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Away from the Ising limit (at nonzero a) the physical 
picture changes qualitatively. First of all, due to the spin- 
flips the spinon is no longer stationary; it may propagate 
through the lattice and has a ±7r/2 momentum in the 
ground state. Second, the spinon-holon interaction V q>q i 
changes. Finally, the holon dispersion is altered as well 
and may acquire some "dressing". The changes in the 
latter, however, should not affect the pairing in any sig- 
nificant way due to the fact that only the holon dispersion 
near the energy minimum matters for it. The holon mass 
renormalization has been analyzed in detail in Ref. [l^ 
and it was found insignificant throughout the anisotropic 
regime < a < 1 . This means that both the "dressing" 
and the holon dispersion changes are minor and should 
not affect pairing. On the other hand, changes in the 
spinon dispersion are qualitative and drastic. At a = 0, 
the spinon may be viewed as a gapped, immobile exci- 
tation with the energy u q = J z /2. With increasing a it 
evolves into a relativistic one, turning completely gap- 
less in the isotropic limit a = 1, where its dispersion is 
= J z (tt/2)\ cosg|. The spinon dispersion for the XXZ 
model at intermediate values of a, shown by solid lines in 
Fig. [21 is known exactly from BA^i see Eq. (Q}. While 
the parameter c/ J z in this equation changes almost lin- 
early between 1/2 and ir/2 as a goes from to 1, the 
parameter k varies from to 1 rather steeply, achieving 
the value of approximately 0.996 at a = 0.5. As a result, 
the spinon gap to s given by lu s = c\f\ — k 2 becomes suf- 
ficiently small already at a ~ 0.5. One can obtain the 
asymptotic behavior for lo s , valid for a > 0.5, and show 
that it approaches zero exponentially in (1 — a)" 1 / 2 as 
a-> 1: 



4c exp 



8(1 - a) 



(15) 



The smallness of the spinon gap may be used to write 
the spinon spectrum in approximate "quasi-relativistic" 
form in this regime: 



(16) 



Another effect of increasing a is a dramatic decrease of 
the spinon's effective mass 



dq 2 



q=Tv/2 



which goes from (4aJ z )~ 1 at a <C 1 to oj s /c 2 at a > 0.5. 
Such a change can be observed in the spectra in Fig. 
[2l where increasing a makes the energy minimum into 
a sharp tip, indicating the mass reduction. Thus, even 
without knowing a specific form of interaction, one can 
anticipate that the spinon-holon binding will be strongly 
affected by such changes in the spinon spectrum. We also 
note that since the spinon becomes much lighter than the 
holon, m s <Cm/t — (2i) _1 , the role of the holon dispersion 
in (fT2")) becomes secondary close to the isotropic limit. 



The remaining question is that of the spinon-holon in- 
teraction. One can analyze the binding problem in the 
small-a limit rigorously. The changes to the holon and 
the AF GSEs are of order 0(a 2 ), while the spinon energy 
changes in the order 0{a): 



Slo„ 



Jz/2 + aj z cos 2q. 



(18) 



One of the consequences of non-zero anisotropy is the 
±7r/2 momentum of the GS of the spinon. This immedi- 
ately implies that the spinon-holon pairing should result 
in a bound state with finite total momentum P = ±7r/2, 
in agreement with the numerical data, shown in Fig. 
[U Since the energy of the system is lowered when the 
AF domain walls associated with the spinon and holon 
pass through each other, the interaction between the two 
can be written as a "contact" attraction of the strength 
V° = — J z /2. Using real-space considerations we find that 
this leads to a direct relation between interaction in the 
momentum space and spinon dispersion. To the order 
0(a) the interaction can be shown to be: 



V q , q - =-u°-(6w q +5u q -)/2. 



(19) 



This equation may be used to derive an analytic expres- 
sion for the binding energy A, exact to the first order in 
a. Substituting (fUD into (JT2|) we get 



Xp{q) 



1 



B 
~2 



where 



A = £xp(</), 

i' 

B = ^2Suj q >Xp(q), 

q' 

E q = E—ep-q—UJq. 



(20) 

(21) 

(22) 
(23) 



After inserting this result for xil) m to (|12[) . dropping 
the higher-order terms in ce, and some algebraic manip- 
ulations, we end up with the following equation for A: 



'-tE 



1 + 2a cos 2q 



(er- 



2a(cos2(7 + 1) ' 



(24) 



(17) Further expansion in a and calculation using the "bare" 



holon energy efe = e 



A which is valid to order 0(a): 

A=A (1-Ca) 
with Ao given by Eq. ([3]) and 

2J 2 



2icosfc, yields an expression for 

(25) 



C = 



^Wt 2 + J 2 (y 16i 2 + J 2 - At 



(26) 



Interestingly, the initial slope of A(a)/A(0) depends only 
weakly on the value of J z : it is bound between C = 4 
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at J z = and C = 2 at J 2 /f > 1 and varies smoothly 
between them. This linear-a result (|25p is shown in Fig. 
[10] with dashed lines. It is in extremely close agreement 
with the numerical data in the small-a regime. 

Having established the form (fT9")l of the spinon-holon 
interaction for small a, we can now try to address the 
question of how might the general form of the interaction, 
valid for any anisotropy a, look like. While there are no 
strict analytical arguments for it, we may formulate a 
number of criteria, which this form must satisfy. First of 
all, the interaction V q ^ q > must be symmetric with respect 
to momenta, V q ^ q > — V q > tq . Second, it must reproduce 
the small a limit (| 19[) as a — > 0. One can also anticipate 
that it should be straightforwardly related to the spinon 
energy, similar to Eq. (fT9"]) . 

Based on these requirements, we propose the following 
form of the interaction in the momentum space: 

Vq,q> = -^JUJqUJq'- (27) 

This is somewhat reminiscent of the electron-phonon 
interaction that is proportional to square-root of the 
phonon energy. Using this Ansatz for V q _ q ', spinon en- 
ergy from BA, and neglecting the changes in the holon 
dispersion (ffc = e°) we arrive at a solution of Eq. (fT2f of 
the form 

X (q) = const x (28) 

leading to the following equation for A: 

l = -Yl^-, \ ■ (29) 

^ A - (ep_ g - e ) - {oj q - uj ) 

Solving this equation numerically yields the complete de- 
pendence of the binding energy A on anisotropy a shown 
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FIG. 13: (Color online). Analytic (At) and numerical (A n ) 
results for the binding energy at small values of a for J z /t = 
4.0. Inset shows the relative difference between the theoretical 
and numerical results. 



as solid lines in Fig. [TO] Not only this equation naturally 
yields our small-a results, but it also provides a very 
close agreement with the numerical data for all values of 
J z and for all a we can access numerically. This pro- 
vides a very convincing a posteriori verification of our 
spinon-holon interaction Ansatz. 

Since we neglect the changes in the holon dispersion, 
the deviation of our theoretical result for A from an ex- 
act answer is expected to occur in order 0(a 2 ). We veri- 
fied that in the small-a limit. Results of the comparison 
of analytic and numerical results arc presented in Fig. 
1131 At a = 0, numerical data obtained using method B 
for the largest system size accessible by ED (19-21 sites) 
agrees with the exact analytical result within the nu- 
merical precision (10 _7 i). However, any small anisotropy 
results in a finite-size effect, linear in a (see inset of Fig. 
1131) . This may be understood in terms of the momen- 
tum space mismatch, discussed earlier: for any finite 
a and finite size L we cannot obtain the "true" GSE 
value for spinon from the numerical simulations, because 
the momentum point corresponding to its lowest energy 
(q = 7r/2) is incommensurate with the available momen- 
tum points. Thus, a finite-size effect of the order 0(1/ L) 
is expected for any finite a. The extrapolated data, on 
the other hand, displays the expected 0(a 2 ) deviation. 

Although we have no formal proof of the validity of 
our interaction Ansatz for all a, the agreement with the 
numerical data makes it very plausible. As the bind- 
ing energy becomes small, it is the long-wavelength fea- 
tures of the dispersions and interaction that determine 
the pairing. One can see from Eq. (|27p that at a — » 1 
the characteristic interaction at low energies is V ~ lu s . 
Thus, within the qualitative picture of pairing in ID, 
both the interaction and the spinon mass become pro- 
portional to the spinon gap to s that tends to zero ex- 
ponentially. One then expects the asymptotic behavior 
A ~ —V 2 m ~ -wj. From Eq. (f2"9")) we can derive such 
an asymptotic expression explicitly: A D(J Z , a)uj^/c 2 , 
where 

f 1, t > J z ; 

V(J z ,a) = \ If 7TC V (30) 
I 2 \4t\n(c/2t)J ' ^ Jz - 

Notably, the exponential behavior of the binding energy 
is determined solely by the asymptotic behavior (|15[) of 
the spinon gap uj s , with the expression in the exponential 
dependent only on a and not on J z /t. Thus, the holon 
energy scale is secondary as it only enters the pref actor. 
Altogether, this explains the quick (exponential) drop-off 

already at intermediate values of a ^ 0.5. 

From this asymptotic expression and Eq. (|29[) one can 
see that the binding energy vanishes in the isotropic limit 
together with the spinon gap. Thus, our spinon-holon in- 
teraction Ansatz also provides a natural and simple ex- 
planation of the non-zero binding at finite q but no bound 
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state at a — 1. This is possible because the interaction of 
the holon with the long-wavelength spinon V q . q > vanishes 
together with the spinon energy. Then, the pairing is not 
strong enough to produce a bound state in the isotropic 
limit. We also find that in the isotropic limit the spinon- 
holon pair wave- function, Eq. (|28|) . is xil) ~ 1/ 1 \f^q- I* 1 
real space, this would correspond to l/y/r spinon-holon 
correlation, exactly the behavior found in Ref. [y. 



IV. CONCLUSIONS 

We have performed extensive analytical and numerical 
studies of an anisotropic version of the t-J model, doped 
with a single hole. Our main result is that the anisotropy 
of the spin-spin interaction leads to an effective attrac- 
tion between the spinon and holon excitations, resulting 
in existence of a spinon-holon bound state. Using the ED 
and DMRG techniques we have numerically estimated 
the binding energy as a function of anisotropy. We have 
described in detail the finite-size effects which arise due 
to various factors and, by examining various ways to mit- 
igate or eliminate them, worked out a procedure for ex- 
trapolation of the finite-size data to the infinite size limit, 
resulting in precise estimates of the binding energy up to 
anisotropy a — 0.5. The resulting numerical values have 
been found to be in excellent agreement with the theory, 
based on Bethe-Salpeter equation. Using the experience 
gained while studying the small anisotropy limit, we have 
formulated the criteria for the form of the spinon-holon 
interaction in momentum space, and proposed a form 
(|27| for it, which results in excellent agreement of ana- 
lytical and numerical results. Finally, we have identified 
the changes in the spinon spectra as the primary factor 
affecting the behavior of the binding energy as a func- 
tion of anisotropy. We have demonstrated that the bind- 
ing energy goes to zero exponentially, as a power of the 
spinon gap, when isotropic limit is approached. This be- 
havior also explains why there is no spinon-holon binding 
in the isotropic t-J model. These results could be tested 
in photoemission experiments in ID spin-chain systems 
with Ising anisotropics. 

We would like to note that the problem we have con- 
sidered is strictly single hole, and its extension to the 
finite-doping case is not trivial. For instance, even in 
the pure Ising limit one could deliberately avoid creating 
any spinons by putting an even number of holes in the 
holon-only states (see Ref. l2lh . However, if spinons are 
present in the system, the interaction between them and 
the holons remains attractive even at a finite hole doping 
and may lead to their binding. 
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APPENDIX A: DIFFERENT TREATMENTS OF 
THE a = PROBLEM 

For pedagogical purposes we describe three different 
analytic ways to determine the bound state energy of the 
t-J z model, doped with a single hole. 

Method 1. The first method we discuss is the exact cal- 
culation of the hole's real-space Green's function. It can 
be accomplished using the expansion in paths ; 22 i 23 or the 
recursion technique 24,25 that is also identical to the Lanc- 
zos method. We use the latter as the most straightfor- 
ward. The bound state energy may be determined as the 
pole of the diagonal element 



1 



H 



(Al) 



where ipi is the state of the system in which the hole is 
located at site i. It may be calculated exactly by noting 
that we may bring the Hamiltonian to a tridiagonal form 
by generating a basis 



where 



\n + l) = H\n)-a n \n)-b 2 n \n-l), 



(n\H\n) 



(n - I n - 1 



(A2) 



(A3) 



It is easy to see that in this basis the diagonal Green's 
function may be represented as a continued fraction: 



GuM = (l 



1 



lo-H 



1 = 



(A4) 



uj — a\ 



In the case of the t-J z model the coefficients in the con- 
tinued fraction have the form 

bi = i 



b\ = 2t 2 

h 2 — h 2 — h 2 — 



and 



01 
(12 



-oj° = -J z /2 

&3 = &4 = ■ • • = 0. 



Thus, we may rewrite (|A4|) as 



where 



E(w) 



oj + oj - 2£(w)' 



uj — £(w) 



(A5) 



(A6) 



(A7) 



(A8) 
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Solving for S(o>) yields the Green's function 

1 



with poles given by 



J z /2T\/w 2 -4t 2 ' 



= ±v/4< 2 + J2/4. 



(A9) 



(A10) 



The energy corresponding to the lowest pole is lower than 
hole's kinetic energy —2t, indicating the presence of a 
bound state. The expression for the binding energy A = 
2t — \Q>\ is equivalent to The quasiparticle residue of 
that pole is given by4 



Z 



J z 
2\Q\ 



Jz 



V 1 ^ 2 + Jz 



(All) 



At large t/J z , in agreement with the naive expectations, 
Z w (J z /2)/2t, the ratio of spinon-holon interaction 
strength to the holon kinetic energy. 

Method 2. Another approach is to consider the immo- 
bile spinon to be an impurity and solve the problem of 
a freely moving hole scattering on it using the T-matrix 
formalism. To that end we write the Hamiltonian of the 
hole as 



-2f y^cos(fc)CfcCfc, 



(A12) 



where c& is the annihilation operator for a hole with mo- 
mentum k. The impurity Hamiltonian may be written 



as 



W = -uJ°clc r , 



(A13) 



i.e. it lowers the energy by cu° = J z /2 if the hole is present 
at the origin site. Transforming it to the momentum 
space after assuming r = yields 



'£4 

k,k> 



Ck>- 



(A14) 



The impurity Hamiltonian is such that each of its matrix 
elements is equal to — u>°. The equation for T-matrix is 



t = h' + h'GqT. 



(A15) 



Inserting complete sets of states and using the fact that 
Go is diagonal and (k\H'\k') = —uj° for arbitrary k, k', 
we get the following equation for the matrix elements: 

(fc|T|fc'> --^ -^°E(p|Go|p)(p|T|fc')- (A16) 



Close examination of this expression reveals that the T- 
matrix is independent of k and k! . Denoting the matrix 
element by T(u>) we obtain 



The integral on the rhs 

S = J2(p\GoH\ P ) 



(A18) 



can be calculated using elementary methods to find 

1 



S = 



Thus 



T(u,) 



Vuj 2 - At 2 
wVw 2 - if 2 



j° + VuJ 2 - it 2 



(A19) 



(A20) 



The exact Green's function of a particle with a scatterer 
is a function of incoming and outgoing momenta k and 
k', given by 



G(k, k') = Go(fc) + G (k)TG Q (k'). 
This yields the following matrix element: 



(A21) 



{k\G(u)\k') 



1 



u> + 2t cos k 



6 k ,k> +T(lu)- 



1 



uj + 2t cos k' 

(A22) 

Transforming to real space by integrating over momenta 
k and k', we get the diagonal Green's function at the 
point of origin: 



G rr (u) = ^(fc|G(w)|fc') 



k.k> 



1 



y/u; 2 - At 2 
1 



At 2 



T( W ) 



j° + Vuj 2 - At 2 



(A23) 



that is equivalent to (|A9[) . This Green's function has the 
poles at the locations given by (|A10|) , so it yields binding 
energy equivalent to ([3]). 

Method 3. Finally, the last approach is based on the 
physical picture of hole decay into a spinon and holon, 
confined to two half-spaces (see Fig. [IJ, and solving the 
Dyson's equation for such a decay exactly. The Dyson 
equation for the hole at origin has the form 



GM 



(GO)-i-SH- 



where 



G°M = -, 



(A24) 



(A25) 



and self-energy E(ct>) may be written in terms of the 
spinon Green's function D(u>) and the holon Green's 

function in a half-space G^ 2 \u>): 



r( W ) = - W ° (l+T( W )5>|GoW|P> . (A17) nuj) = 2t 2 [ d ^D{u')G^\u-u'). (A26) 

\ v J J 2n 
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Since 



D(u) 



(A27) 



the integral is readily done, leading to the expression for 
self-energy 



(A28) 



The Green's function (|A24|) can then be written in terms 
of the shifted frequency u> = w — w as 



G(u>) 



UJ + UJ 



E(S) 



(A29) 



To calculate G^ 1 (£>), that is the Green's function of a 
free hole subject to a "hard- wall" boundary condition at 
the origin, we note that it can be expressed as the anti- 
symmetric part of the hole's Green's function Gh{oj) in 
the entire space. For example, in the coordinate repre- 
sentation we obtain: 

g£ /2) (x,x';uj) = G h (x,x';w) - G h (x,-x';u). (A30) 



Here Gh{x,x' ;uj) is just an inverse Fourier transform of 
the Green's function of a free hole: 



Gh{x,x';u) 



dp e 



i(x—x')p 



2ir ui + 2t cosp 



(A31) 



It depends only on the difference of the coordinates x—x', 
as expected for a translationally invariant system. We are 
interested in the diagonal matrix element of (|A30[) : 

Gi 1/2) (x; W ) = G^ /a \x,x;u) = 

= G h (x-x;u)-G h (x + x;u;).(A32) 



Using (TA3T|l it may be readily found to be 

2G^ 1/2) (u)=w± Vlo 2 - At 2 , (A33) 

again leading to the expression for G(u>) equivalent to 
Eqs. {El and (|A23| and expression for A, equivalent to 
©■ 
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